# =====================================================
# Replication Material: Main Results
# König et al. (2019)
# Comparative Politics and Causal Evaluation of Structural Reforms: The Case of the UK National Minimum Wage Introduction
# =====================================================

# start log file
sink("./log_main.txt")

rm(list = ls())
# Install packages needed for analysis
p_needed <- c("devtools", "ggplot2", "dplyr", "tidyr", "xtable") 
packages <- rownames(installed.packages())
p_to_install <- p_needed[!(p_needed %in% packages)]
if (length(p_to_install) > 0) {
  install.packages(p_to_install)
}

# Main analysis with gsynth is based on version 1.0.2
library(devtools)
install_version("gsynth", version = "1.0.2", repos = "http://cran.us.r-project.org")

# load packages
library(gsynth)
library(ggplot2)
library(dplyr)
library(tidyr)
library(xtable)

# =================
# load data
# =================
rm(list=ls())
load("01_data/sample.RData")

# =================
# load plotting function
# =================
# contains minor adjustments to the plotting function in gsynth to make it look prettier
source("function_plot_gsynth.R")

# =================
# produce results
# =================

# Figure 1
# ==========
# run main analysis
out.1 <- gsynth(yurate ~ mw 
                         , data = sample, 
                         index = c("ccode","Year"), force = "two-way",
                         CV = TRUE, r = c(0, 5), se = TRUE, 
                         inference = "parametric", nboots = 1000,
                         parallel = F, EM=TRUE, seed = 1021)

# Counterfactual plot
plot(out.1, type = "counterfactual", xlab="", ylab = "Youth Unemployment Rate (%)", main = "", raw = "none") + 
  theme_classic(20, base_family = "serif") + geom_vline(xintercept = 1998,linetype="dashed") +  
  theme(legend.position="bottom") + geom_vline(xintercept = 2004,linetype="dotted") +
  geom_vline(xintercept = 2008,linetype="dotted") +
  scale_x_continuous(breaks=c(1985, 1990,1995,2000,2005,2010)
                     , labels = c("1985","1990","1995","2000","2005","2010")
  ) + ylim(0,25)

ggsave("02_figures_tables/Figure1.pdf"
       , width=(6.3228344444445)*1.5, 
       height=6.3228344444445, units="in")

# Figure 2a
# ==========

# Gap plot
plot(out.1, type="gap", xlab = "", ylab = "Gap Youth Unemployment (% points)",  main = "") + 
  theme_classic(20, base_family = "serif") + geom_vline(xintercept = 1998, linetype="dashed") +  
  geom_hline(yintercept = 0,linetype="dashed") + ylim(-12,20) +
  scale_x_continuous(breaks=c(1985, 1990,1995,2000,2005,2010)
                     , labels = c("1985","1990","1995","2000","2005","2010")
  )
ggsave("02_figures_tables/Figure2a.pdf"
       , width=(6.3228344444445)*1.5, 
       height=6.3228344444445, units="in")

# Figure 2b
# ==========

# run analysis with Eastern Immigration as share of active labor force as control variable
out.2 <- gsynth(yurate ~ mw + log(east.labforce)
                           , data = sample, 
                           index = c("ccode","Year"), force = "two-way",
                           CV = TRUE, r = c(0, 5), se = TRUE, 
                           inference = "parametric", nboots = 1000,
                           parallel = F, EM=TRUE,seed = 1001)


# Gap plot
plot(out.2, type="gap", xlab = "", ylab = "Gap Youth Unemployment (% points)",  main = "") + 
  theme_classic(20, base_family = "serif") + geom_vline(xintercept = 1998,linetype="dashed") +  geom_hline(yintercept = 0,linetype="dashed")+ ylim(-12,20) +
  scale_x_continuous(breaks=c(1985, 1990,1995,2000,2005,2010)
                     , labels = c("1985","1990","1995","2000","2005","2010")
  )
ggsave("02_figures_tables/Figure2b.pdf"
       , width=(6.3228344444445)*1.5, 
       height=6.3228344444445, units="in")

# Figure 2c
# ==========

# run analysis with delta GPD as control variable
out.3 <- gsynth(yurate ~ mw + DeltaGDP
                           , data = sample, 
                           index = c("ccode","Year"), force = "two-way",
                           CV = TRUE, r = c(0, 5), se = TRUE, 
                           inference = "parametric", nboots = 1000,
                           parallel = F, EM=TRUE,seed = 1008)

# Gap plot
plot(out.3, type="gap", xlab = "", ylab = "Gap Youth Unemployment (% points)",  main = "") +
  theme_classic(20, base_family = "serif") + geom_vline(xintercept = 1998,linetype="dashed") +  geom_hline(yintercept = 0,linetype="dashed") + ylim(-12,20) +
  scale_x_continuous(breaks=c(1985, 1990,1995,2000,2005,2010)
                     , labels = c("1985","1990","1995","2000","2005","2010")
  )
ggsave("02_figures_tables/Figure2c.pdf"
       , width=(6.3228344444445)*1.5, 
       height=6.3228344444445, units="in")

# Figure 2d
# ==========

# run analysis with Eastern Immigration as share of active labor force and delta GDP as control variables
out.4 <- gsynth(yurate ~ mw + log(east.labforce) + DeltaGDP
                           , data = sample, 
                           index = c("ccode","Year"), force = "two-way",
                           CV = TRUE, r = c(0, 5), se = TRUE, 
                           inference = "parametric", nboots = 1000,
                           parallel = F, EM=TRUE,seed = 1136)


# Gap plot
plot(out.4, type="gap", xlab = "", ylab = "Gap Youth Unemployment (% points)",  main = "") + 
  theme_classic(20, base_family = "serif") + geom_vline(xintercept = 1998,linetype="dashed") +  geom_hline(yintercept = 0,linetype="dashed")+ ylim(-12,20) +
  scale_x_continuous(breaks=c(1985, 1990,1995,2000,2005,2010)
                     , labels = c("1985","1990","1995","2000","2005","2010")
  )

ggsave("02_figures_tables/Figure2d.pdf"
       , width=(6.3228344444445)*1.5, 
       height=6.3228344444445, units="in")

# Table 1
# ==========
tab.reg <- rbind(
  # ATT
  cbind("Average Treatment Effect NMW", round(out.1$att.avg,2), round(out.2$att.avg,2), round(out.3$att.avg,2), round(out.4$att.avg,2))
  # se
  ,cbind("", paste("(",round(out.1$est.avg[2],2),")", sep = ""), paste("(",round(out.2$est.avg[2],2),")", sep = ""), paste("(",round(out.3$est.avg[2],2),")", sep = ""), paste("(",round(out.4$est.avg[2],2),")", sep = ""))
  # EE
  ,cbind("EE immigration (\\% of active labor force)", NA, round(out.2$est.beta[1],2), NA, round(out.4$est.beta[1,1],2))
  # se
  ,cbind("" , NA, paste("(",round(out.2$est.beta[2],2),")", sep = ""), NA, paste("(",round(out.4$est.beta[1,2],2),")", sep = ""))
  # DeltaGDP
  ,cbind("$\\Delta$ GDP pc", NA, NA, round(out.3$est.beta[1],2), round(out.4$est.beta[2,1],2))
  # se
  ,cbind("", NA, NA, paste("(",round(out.3$est.beta[2],2),")", sep = ""), paste("(",round(out.4$est.beta[2,2],2),")", sep = ""))
  # State FE
  ,cbind("State FE", t(rep("x",4)))
  # Year FE
  ,cbind("Year FE", t(rep("x",4)))
  # Unobserved Factors
  ,cbind("Unobserved Factors",out.1$r.cv, out.2$r.cv, out.3$r.cv, out.4$r.cv)
  # N
  ,cbind("N",dim(out.1$Y.dat)[1]*dim(out.1$Y.dat)[2]
         ,dim(out.2$Y.dat)[1]*dim(out.2$Y.dat)[2]
         ,dim(out.3$Y.dat)[1]*dim(out.3$Y.dat)[2]
         ,dim(out.4$Y.dat)[1]*dim(out.4$Y.dat)[2]
  )
  # Control States
  ,cbind("Control States",out.1$Nco, out.2$Nco, out.3$Nco, out.4$Nco)
)

# save output
print(xtable(tab.reg
             , digits = c(0, 0, 2, 2, 2, 2)
             , caption = "The effect of the UK NMW on youth unemployment rate 1999-2012"
             , label = "tab:results"
)
, include.rownames = FALSE
, include.colnames = FALSE
, caption.placement = "top"
, sanitize.text.function = force
, table.placement = "t!"
#, hline.after = c(-1,0,nrow(x))
, add.to.row = list(pos = list(0), command = c("&   (1) & (2) & (3) & (4) \\\\"))
, file = "02_figures_tables/Table_1.tex"
)

# ======================
# ======================
# close log file
sink()
